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Abstract 

The Helmholtz equation arises in many applications, such as seismic 
and medical imaging. These application are characterized by the need to 
propagate many wavelengths through an inhomogeneous medium. The 
typical size of the problems in 3D applications precludes the use of direct 
factorization to solve the equation and hence iterative methods are used 
in practice. For higher wavenumbers, the system becomes increasingly 
indefinite and thus good preconditioners need to be constructed. In this 
note we consider an accelerated Kazcmarz method (CGMN) and present 
an expression for the resulting iteration matrix. This iteration matrix can 
be used to analyze the convergence of the CGMN method. In particular, 
we present a Fourier analysis for the method applied to the ID Helmholtz 
equation. This analysis suggests an optimal choice of the relaxation pa- 
rameter. Finally, we present some numerical experiments. 



1 Introduction 

The Helmholtz equation arises in many applications such as seismic and medical 
imaging. These application are characterized by the need to propagate many 
wavelengths through an inhomogeneous medium. The scale of the problems in 
3D application precludes the use of direct factorization to solve the equation 
and hence iterative methods are used in practice. For higher wavenumbers, the 
system becomes increasingly indefinite and thus good preconditioners need to be 
constructed. For an overview of preconditioning techniques for the Helmholtz 
equation we refer to [5] and references cited therein. 

In this note, we consider an accelerated Kazcmarz method (CGMN, pQ) for 
solving the Helmholtz equation [4]. We present an expression of the iteration 
matrix that allows us to analyze its convergence behaviour. 

We proceeds as follows. First, we give a brief overview of the CGMN method 
and its relation to SSOR. This gives us an expression for the iteration matrix 
in terms of the original matrix. Then, we present the Fourier analysis for a ID 
finite-difference discretization of the Helmholtz equation and we present some 
numerical examples. Finally, we draw conclusions and discuss future work. 



1.1 The CGMN method 



The CGMN method pQ relies on a symmetric Kaczmarz (SKACZ) sweep in 
conjunction with the conjugate gradient (CG) method. 

The Kaczmarz method solves a system of TV equations, Au — s, by cyclically 
projecting the iterate onto rows of the matrix [5] 

u:=u + uj(s t -a l T u)a l /\\a. l \\l, i = l...N, (1) 

where denotes the i-th row of A as column vector and < to < 2 is a 
relaxation parameter. Note that we may also let u> vary per row. Introducing 
the matrices Q b = [I — oj&i&i T /\ \&i\ \%) , we may write this iteration as 

u := QjU + waj/Hajll^ (2) 

A double sweep through the matrix (from row 1 to N and back) can then be 
denoted by 

u:=Qu + i?s, (3) 

where Q = QoQi ■ ■ . Qn-iQn-i ■ ■ ■ Qo and R contains all the factors multiply- 
ing b. It is easily verified that the matrices Qi are symmetric and have rank 
1 with eigenvalue 1 — uj. It follows that Q is symmetric and has eigenvalues 
€ [—1, 1]. Hence, (/ — Q) is symmetric and positive semi-definite and we can 
use CG to solve the equivalent system 

(I-Q)u = Rs. (4) 

The question is, how does this new system behave in terms of w? In the next 
section, we will first derive a convenient expression for the matrix Q in terms 
of the matrix A, and in the subsequent section we present a Fourier analysis of 
the error propagation applied to the ID Hclmholtz equation. 



1.2 Relation to SSOR 

It is well-known that a SKACZ sweep on the system ^4u = s is equivalent to a 
SSOR sweep on the normal equations AA T x = s, u = A T x pQ. Splitting the 
matrix AA T into its diagonal, strictly lower and upper triangular parts D, L 
and L T such that AA T = D + L + L T , the SSOR iteration can be represented 
as follows [7 

x:=Gx + Hs, (5) 

where 

G = {D + ujL t )- 1 {{1-uj)D-loL){D + ujL)- 1 {{1-uj)D-loL t ), (6) 
H = oj(2-oj)(D + ojL t )- 1 D(D + ujL)- 1 . (7) 
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Since the SKACZ and SSOR iterations are equivalent, the matrices Q, R and 
G, H are related via 

QA T = A T G, (8) 
R = A T H. (9) 

Next, we note that we may write G in terms of H and A as 

G = I — HAA T , (10) 

from which it follows that we can write 

Q = I — A T HA. (11) 

2 Fourier analysis of the error propagation 

We study the properties of the preconditioner by looking at the error propagation 
matrix. This matrix relates errors in u when doing a standard Richardson 
iteration on the preconditioned system: 

u k+1 = (I - (I - Q))u k+1 + Rb. (12) 

In terms of the error e^+i = — u*, we get 

efe+i = Qe k , (13) 

Similarly, the residual = Ae k is propagated by 

r fe+1 = G T r k . (14) 

We analyze the behaviour of the SKACZ preconditioner through a Fourier 
analysis of the error propagation matrix [21 [6] . We decompose the error into its 
Fourier modes ej(6) — e^ 9 and want to find the amplitude a(8) such that 

Qe(6) = a(6)e(6), (15) 

Using eq. Q and ( 11 ) we can factorize the amplitude function in terms of those 
of the matrices A (oi), (D + ujL) (02), D (03) and (D + u>L T ) (04): 

o = l-w(2-w)^. (16) 
a 2 a 4 

Note that the matrix G has the same amplitude function. 
2.1 ID case 

We consider a simple finite difference discretization of the the ID Helmholtz 
equation (k 2 + d 2 /dx 2 )u = son the domain ft = [0, 1] with Dirichlet boundary 
conditions. The resulting system of equations is denoted by 

Au = s, (17) 
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where 

an = k 2 -2h~ 2 , (18) 

<H,i±i = h~ 2 , (19) 

and h = 1/(JV+1) denotes the gridspacing and N denotes the number of interior 
gridpoints. The matrix A is symmetric and has eigenvectors 

v" = sm(mrih), n = l,2,...,N (20) 

with corresponding eigenvalues 

A„ = k 2 + 2h- 2 (cos(nirh) - 1). (21) 

The entries of the matrices D and L are given by da — | 1 1| = 7 2 + 2h~ A , 
li t i-i = af ai_i = 2h~ 2 ^ and £j,j_2 = /i~ 4 where 7 = fc 2 — 2/i~ 2 . 
It is now readily verified that 

«i = 7 + 2/i~ 2 cos6», (22) 

a 2 = 7 2 + 2/i- 4 + w/i- 2 (27exp-^+/i- 2 exp- 2,,? ), (23) 

a 3 = 7 2 + 2^ 4 , (24) 

a 4 = 7 2 + 2^ 4 + cj/i- 2 (27exp +je +/i- 2 exp +2i9 ). (25) 



Substituting these expressions in eq. ([16]) yields 
a(9, ui, h) = 1— 

w(2- w)^( 7 + : 



/3 2 + 2/3w(27/i- 2 cos (9 + cos 2(9) + w 2 /i- 4 (4 7 2 + 4 7 /i~ 2 cos + h^) ' 

(26) 

where /3 = 7 2 + 2/i~ 4 . 

To ensure a certain amount of gridpoints per wavelength, n g , we let kh = 
2tt /rig. Note that in this case, the amplitude only depends on n g and oj. The 

a) shows the 
(b) we show 

maxa(i-a(g)) ^ w j 1 j c i 1 j s r0 ughly equivalent to plotting the condition number of 



question is, what is the optimal u for a given n g . Figure 
amplitude as a function of and ui for n g = 10. In Figure 



miny (1 — a{0)) ' 

I — Q, and this suggests an optimal value of oj ss 1.5 for n g = 10. Figure [4] (c) 
shows the optimal w as a function of n g . We can use this curve to adapt uj to 
the (local) wavenumber. 



3 Numerical experiments 

We present some numerical experiments that verify the results discussed above. 
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3.1 ID 



In the first experiment we solve the Hclmholtz equation for various wavenumbers 
k, with a fixed number of n g = 10 gridpoints per wavelength. FigureBl(a) shows 
the number of CG iterations needed to converge to a tolerance of 10~^or various 
ui. The predicted optimal value (cf. figure [4] (c)) is also indicated and nicely 
coincides with the optimal ui in terms of the iteration count. 

In the second experiment, we use the same values of k, but choose the 
gridspacing based on the highest wavenumber used. This way, the number of 
gridpoints per wavelength is much higher for the lower wavenumbers. The result 
is shown in Figure [4] (b) . Again, the predicted optimal value of u is close to the 
empirical optimum. 



3.2 2D 

We use a medium that consist of a high-wavenumber anomaly embedded in 



r e [4| 



(a). The 
(b)) is de- 



a background with constant wavenumber ko, as depicted in figure 
wavefield resulting from an incident planewave exp(ifcox) (figure I 
picted in figure |4](c). The convergence histories using either a constant oj = 1.5 
(chosen according to the highest wavenumber) or uj — 1.95 (chosen according 
to the lowest wavenumber) or a spatially varying uj (corresponding to the local 
wavenumber) are shown in figure [4] (c) . It appears that for a medium with a 
large contrast it may indeed be beneficial to vary the relaxation parameter spa- 
tially. Interestingly enough, the convergence for uj = 1.95 is much faster in the 
beginning, suggesting perhaps that it may be beneficial to change the relaxation 
parameter per iteration. 



4 Conclusion and future work 

We have presented a Fourier analysis of the CGMN method applied to the 
ID Hclmholtz equation. The amplitude function of the error propagation ma- 
trix suggests that the optimal relation paramater depends on the number of 
gridpoints per wavelength n g = 2ir/{kh). This observation was confirmed by 
counting the number of CGMN iterations needed to actually solve the system 
for various values of uj and n g . The same optimal relation can be used to adapt 
the relaxation parameter to the local wavenumber in case of inhomogeneous me- 
dia. A numerical example on a 2D medium with a very high contrast suggests 
that this might indeed improve convergence. 

Future research is aimed at extending this analysis to 2D/3D Helmholtz 
equations and varying media, as well as the parallel extension of the CGMN 
method presented by [3] . An optimal strategy to adapt the relaxation parameter 
to the medium and perhaps vary it per iteration is of particular interest. 



5 



Acknowledgments 



The author thanks Dan and Rachel Gordon for valuable discussions on this 
topic. This work was in part financially supported by the Natural Sciences and 
Engineering Research Council of Canada Discovery Grant (22R81254) and the 
Collaborative Research and Development Grant DNOISE II (375142-08). This 
research was carried with support from the sponsors of the SINBAD consortium. 

References 

[1] Ake Bjorck and Tommy Elfving. Accelerated projection methods for comput- 
ing pscudoinverse solutions of systems of linear equations. BIT, 19(2):145- 
163, June 1979. 

[2] A. Brandt. Multi-level adaptive solutions to boundary-value problems. 
Mathematics of Computation, 31:333-390, 1977. 

[3] O. G. Ernst and M. J. Gander. Why it is Difficult to Solve Helmholtz Prob- 
lems with Classical Iterative Methods, volume 83 of Numerical Analysis of 
Multiscale Problems. 2011. 

[4] Dan Gordon and Rachel Gordon. Robust and highly scalable parallel solu- 
tion of the Helmholtz equation with large wave numbers. Journal of Com- 
putational and Applied Mathematics, 237(1):182-196, January 2013. 

[5] S. Kaczmarz. Angenaherte Auflosung von Systemen linearer Glcichungen. 
Bulletin International de 1' Academic Polonaise des Sciences et des Lettres, 
35:355-357, 1937. 

[6] Rob Kettler. Analysis and comparison of relaxation schemes in robust 
multigrid and preconditioned conjugate gradient methods. In W. Hack- 
busch and U. Trottenberg, editors, Multigrid Methods, volume 960 of Lecture 
Notes in Mathematics, pages 502-534. Springer Berlin / Heidelberg, 1982. 
10.1007/BFb0069941. 

[7] Yousef Saad. Iterative methods for sparse linear systems. SIAM, 2 edition, 
2003. 



6 
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(b) 




25 30 



Figure 1: (a) Amplitude (cf. eq. 26 1 as a function of 9 and to for a fixed number 
of gridpoints per wavelength n g — 10. (b) Shows the 'condition number' as a 
function of uj for n g = 10. The optimal ui ss 1.5. (c) shows the optimal w as a 
function of n„. 




wavenumber [1/m] wavenumber [1/m] 



(a) (b) 

Figure 2: Number of CG iterations needed to converge to a tolerance of 10~ 6 for 
various co with either a fixed number of gridpoints per wavelength n g = 10 (a) 
or a fixed gridspacing based on the highest wavenumber used (b). The predicted 
optimal to is indicated by a dashed line and coincides with the lowest iteration 
count. 
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Figure 3: (a) wavenumber profile, (b) incident plane wave, (c) resulting wavefield 
and (d) convergence histories for u) = 1.5 (blue), ui — 1 .95 (green) and a spatially 
varying cj (red). 
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